An Upwind-Biased Space 
Marching Algorithm for 
Supersonic Viscous Flow 


Francis A. Greene 





NASA 

Technical 

Paper 

3068 


1991 


An Upwind-Biased Space 
Marching Algorithm for 
Supersonic Viscous Flow 


Francis A. Greene 
Langley Research Center 
Hampton, Virginia 


fWNSA 

National Aeronautics and 
Space Administration 

Office of Management 

Scientific and Technical 
Information Division 



Symbols 

A 

a, b 

Cn 

c 

D 

E 

E' 

E 

e 

e' 

F 

G 

H 

H 

h 

i,j,k 

L 

I 

lx-, lyi lz 
M 
Af f 

M 
M 1 

to 

m x ,m y ,m z 

n 

Tl y-> Tiz 

P 

P 

P f 

Pr 

Q 

q f 

Qx 

Qy 


inviscid flux Jacobian matrix 

inviscid flux volume weighting parameters 

Courant number 

nondimensional sound speed 

flux difference 

nondimensional total energy per unit volume 
dimensional total energy per unit volume, kg/m-s 2 
flux vector in £ direction 

nondimensional internal energy per unit mass 
dimensional internal energy per unit mass, m 2 /s 2 
flux vector in 7} direction 
flux vector in £ direction 
nondimensional total enthalpy per unit mass 
flux vector normal to cell face 

Cartesian flux tensor 

indices of cell centers in £, r], C directions 

normalizing reference length, m 

unit vector tangent to cell face 

components of l in x, t/, z directions 

Mach number 

streamwise Mach number 

right inviscid Jacobian eigenvector matrix 

left inviscid Jacobian eigenvector matrix 

unit vector tangent to cell face 

components of m in x, y, 2 directions 

unit vector normal to cell face 

components of n in x, y, z directions 

dummy index for i, or k 

nondimensional pressure 

dinensional pressure, N/m 2 

Prandtl number 

heating rate 

dimensional heating rate, w/m 2 
nondimensional heating rate in x direction 
nondimensional heating rate in y direction 

iii 




PRECEDING PAGE BLANK NOT FILMED 


Qz 

q 

Re 

R n 

S' 

s 

T 

T' 

t 

U 

U 

VL 

V 


W 

U) 

u/ 

x,y,z 

& 

7 

A 

b 

A 

P 




P 

a 


n 

Subscripts: 

/ 


nondimensional heating rate in ^ direction 

dependent variable vector 

Reynolds number 

nose radius 

velocity, m/s 


nondimensional Cartesian velocity vector 
nondimensional temperature 
dimensional temperature, K 
time, s 

nondimensional velocity in direction of n 

nondimensional Cartesian velocity component in x direction 

dimensional Cartesian velocity component in x direction, m/s 

nondimensional velocity in direction of l 

nondimensional Cartesian velocity component in y direction 

dimensional Cartesian velocity component in y direction, m/s 

nondimensional velocity in direction of fa 

nondimensional Cartesian velocity component in w direction 

dimensional Cartesian velocity component in z direction, m/s 

nondimensional Cartesian coordinates 


dimensional Cartesian coordinates, m 
= 7-1 

ratio of specific heats 
spatial change in quantity, ( )/? — ( )i 
time change in quantity, ( ) n+1 — ( ) n 
inviscid Jacobian eigenvalue matrix 
nondimensional viscosity 
viscosity, kg/m-s 


2 



computational coordinates 
nondimensional density 
dimensional density, kg/m 3 
nondimensional cell face area 
shear stress 

nondimensional cell volume 


inviscid quantity 


iv 



L 

R 

V 

oc 

Superscripts: 
n, * 

T 

+ 


quantity at cell center to left of cell face 
quantity at cell center to right of cell face 
viscous quantity 
free-stream value 

time level 
transpose 

quantity associated with positive eigenvalue 
quantity associated with negative eigenvalue 


v 



Summary 

The Langley Aerothermodynamic Upwind Relaxation Algorithm, which globally iterates 
the thin-layer Navier-Stokes equations, is converted to an implicit space marching algorithm. 
The original code is used to generate the starting solution for its newly developed counterpart 
Because the two algorithms differ only in the computation of the streamwise flux, an added 
degree of consistency between the starting and marching codes is achieved. The space marching 
algorithm is first- or second-order accurate with Roe’s upwind differencing or symmetric total 
variation diminishing differencing, respectively. The upwind formulation of the governing 
equations inherently limits the pressure difference across a cell and avoids the need tor a 
limiting variable placed directly within the governing equations. Each cross-flow plane is locally 
iterated in pseudo time until converged, then marched in space to the next station where the 
convergence process is repeated. This procedure gives a complete solution in a single sweep 
over the geometry. The algorithm is tested on a sphere-cone geometry at 0 incidence and a 
angle of attack on a geometry which models the windward surface of the Space Shuttle whiter 
Center-line surface pressure comparisons for the sphere-cone are made between the globa an 
space marching algorithms. In the overexpansion recompression region of the sphere-cone, 
the space marching algorithm was very sensitive to the step size. The effects of streamwise 
grid refinement on solution quality in the overexpansion recompression region of the sphere- 
cone are discussed. In addition, computational results for surface heating are compared with 
ground-based experimental data for both the sphere-cone and the shuttle-like geometries. 

1. Introduction 

Upwind methods are well suited for the numerical simulation of hypersonic flows, since 
their robust nature enables them to capture strong bow shocks without oscillations. Initially, 
upwind techniques such as the flux-difference split methods developed by Godunov (ret. 1) 
Osher (ref. 2), and Roe (ref. 3), and the flux- vector split methods attributed to Steger and 
Warming (ref. 4) and Van Leer (ref. 5) were developed for the Euler equations. These methods 
have since been used with the Navier-Stokes equations (refs. 6 to 9). More recently, the flux- 
difference split method of Roe has been used with space marching techniques to combine the 
efficiency of space marching with the shock capturing capability of upwind methods in a single 

code (refs. 10 and 11). A 

Prior to upwind differenced algorithms, space marching codes that solved the parabolized 

Navier-Stokes (PNS) equations used solution algorithms based on central difference methods 
(refs. 12 and 13). Unlike solution algorithms based on upwind differencing, centrally differenced 
algorithms can introduce spurious oscillations in values across shock waves. Because of these 
oscillations, centrally differenced algorithms usually have dissipation added into the solution 
procedure to maintain stability. Upwind methods are naturally dissipative and as such require 

no additional dissipation. , 

Even though flow conditions are more restrictive for space marching algorithms compared 
with globally iterated algorithms, the space marching aspect brings a higher degree of efficiency 
when compared with global algorithms. Efficiency in terms of computer time is enhanced 
with space marching techniques. This results in quicker turnaround times and enables space 
marching algorithms to be used on a wider range of computers. 

The thin-layer Navier-Stokes (TLNS) equations are a viable and practical alternative to 
the full Navier-Stokes equations when streamwise viscous effects are negligible. With the 
same viscous terms as the PNS equations, the TLNS equations retain the unsteady terms. 
Even though the TLNS equations retain time-dependent terms, they can be spatially marched 
in the spirit of the PNS equations as long as the appropriate provisions are made at the 
outflow boundary. The globally iterated thin-layer code provides the initial data needed tor 
the marching algorithm. Data at the station where the solution is to be obtained are locally 
iterated in pseudo time until they are converged. Once converged, the solution is marched 
to the next station where the local iteration process is repeated. Unlike a globally iterated 
solution, once a station has been converged one time it is never recomputed. This solution 
procedure has been demonstrated with flux-vector splitting for mviscid flow by Walters and 


Dwoyer (ref. 14) and for viscous flow by Newsome, Walters, and Thomas (ref. 15). In both 
cases, a globally iterated algorithm was modified to space march. Walters and Dwoyer handle 
the subsonic region with a hybrid procedure. Global iterations are first performed identifying 
the subsonic regions, and the solution is then globally iterated in the subsonic regions and space 
marched in the supersonic regions. Newsome et al. use the approach of Vigneron, Rakich, and 
lannehill (ref. 16) to avoid departure solutions. 

The Langley Aerothermodynamic Upwind Relaxation Algorithm (LAURA) (ref. 6) has 
been converted to a space marching algorithm. This report details the modifications made 
to the original algorithm and presents comparisons with experimental data. In converting the 
a gorithm to perform space marching, the original formulation was not disturbed so that in the 
future a single formulation could be used to calculate globally iterated solutions, space-marched 
solutions, or a hybird of both. 

To accurately describe the physics of the flow, the flux (viscous and inviscid) at a cell 
face must be formulated to represent the proper zone of dependence. Second-order accurate 
central difference approximations to the dissipative terms in the viscous flux vector can be 
used to adequately describe the viscous zone of dependence but cannot be used to approximate 
t e terms in the inviscid flux vector such that the inviscid zone of dependence is properly 
reflected. To reflect the inviscid zone of dependence, the motion of waves relative to a cell 
face must be properly accounted for. The propagation of information within the inviscid zone 
of dependence, in terms of speed and direction, is reflected in the eigenvalues of the unsteady 
inviscid flux vector Jacobian. The unsteady Jacobian is considered since the equations to be 
solved will be iterated in time. The space marching algorithm employs the flux-difference 
split method of Roe to define an upwind-biased inviscid flux at cell face in the streamwise as 
well as the cross-flow directions. Inherent in an upwind-biased flux is the knowledge that it 
contains the amount of upstream and downstream information needed to correctly describe the 
inviscid zone of dependence. As a consequence, the flow physics is correctly modeled in one 
dimension and approximately modeled in two and three dimensions. Supersonic inviscid fluxes 
with eigenvalues which are all positive are comprised solely of upstream information, whereas 
fluxes with mixed-sign eigenvalues contain both upstream and downstream information. 

The streamwise pressure gradient propagates information upstream through the subsonic 
portion of the boundary layer. If the elliptic nature of this upstream movement of information, 
which is associated with the negative eigenvalue, is not properly treated, departure solutions will 
result (ref. 17). To avoid departure solutions, the governing equations can be globally iterated 
where the sweep direction is alternated. This approach is typically used when time-dependent 
P resent but can also be used on equation sets without time-dependent terms such as 
the PNS equations (refs. 18 to 20). By performing multiple relaxation sweeps over the body, 
information is allowed to propagate through the subsonic portion of the boundary layer without 
causing the solution to diverge. A consequence of this method is the increased computational 
time to perform multiple sweeps compared with a single sweep method. 

To prevent single sweep space marching methods from being ill-posed, the aforementioned 
elliptic effects must be suppressed. The suppression requires that the influence of the negative 
eigenvalue be limited, this is usually accomplished through a modification of the streamwise 
pressure gradient. To use single sweep space marching methods, the streamwise pressure 
gradient must be favorable (negative) or not so adverse as to cause reversed flow, and the 
streamwise flow should be supersonic everywhere except in the boundary layer. Through an 
eigenvalue analysis of the PNS equations, Vigneron et al. (ref. 16) determined that a fraction of 
the streamwise pressure gradient could be directly included in the PNS equations and the 
remaining portion would either have to be neglected or explicitly imposed. The fraction 
included directly in the PNS equations represents the amount of the pressure gradient associated 
with the positive eigenvalues, whereas the remaining portion is associated with the negative 
eigenvalues. Another approach incorporated within single sweep space marching PNS solvers 
to handle departure solutions is attributed to Schiff and Steger (ref. 21). Their method imposes 
the pressure gradient from the first supersonic point outside the boundary layer throughout 
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the subsonic boundary layer. Lubard and Helliwell (ref. 22) tried lagging the gradient and 
calculating it explicitly. This worked with moderate success. r 

In the present report, the full pressure in the momentum flux is retained. When the 
upwind-biased flux is computed, the pressure in the momentum flux consists of upstream and 
downstream information. This flux, when used in conjunction with a proper specification for 
a supersonic outflow boundary condition, results in the streamwise pressure difference across a 
cell being a fraction of the total change in pressure across the cell. This approach is similar to 
Vigneron’s; however it does not require an explicitly defined limit on the streamwise pressure 
gradient in the governing equations. The limit on the pressure gradient used in Vigneron’s and 
in this approach is different. The differences and similarities are described in this report. 

The space marching algorithm is applied on a 15° sphere-cone geometry at zero incidence 
and at angle of attack on a geometry which models the windward surface of the Space Shuttle 
orbiter (ref. 23). This geometry is referred to as the HALIS shuttle geometry. The sphere- 
cone solutions for surface heating are compared with ground-based experimental data, whereas 
surface pressures are compared with pressures from the globally iterated algorithm. Streamwise 
grid refinement in the overexpansion recompression region of the cone is performed to determine 
the space marching algorithm's dependence on step size in regions of adverse pressure gradients. 
To provide a rigorous test for the new algorithm, space marching solutions on the geometry 
which models the windward surface of the Space Shuttle were computed. Computed surface 
heating distributions at an angle of attack of 10° are compared with ground-based experimental 

data. 


2. Algorithm Development 

The thin-layer Navier-Stokes equations for a compressible, perfect gas are employed in the 
present work. The thin-layer equations are derived from the full Navier-Stokes equations by 
neglecting derivatives in the streamwise and circumferential directions found in the viscous 
terms. This equation set provides adequate solutions for flows in which viscous gradients in 
the streamwise and circumferential directions are small compared with viscous gradients in the 
normal direction. Under these conditions, expending computer resources to resolve viscous 
effects parallel to the body is an inefficient use of computer memory and time, since their effect 
on the converged solution is insignificant. 

The governing equations upon which the space marching algorithm was based were taken 
from reference 6. The program mentioned in reference 6 is the Langley Aerothermodynamic 
Upwind Relaxation Algorithm (LAURA). The inviscid fluxes are either first-order accurate with 
an upwind technique attributed to Roe (ref. 3) or second-order accurate with the symmetric 
total variation diminishing (STVD) scheme of Yee (ref. 24). Possible entropy violations are 
alleviated by a procedure of imposing a minimum on the absolute value of the eigenvalues ol 
the inviscid Jacobian matrix. This fix was suggested by Harten (ref. 25). The terms in the 
viscous flux vector are approximated with second-order central differences. 

The integral form of the governing equations can be expressed as 

J J J qf dQ + j J ti n da = 0 ( 2 -l) 


where q, is the time rate of change of the dependent variable vector at a cell center, h is the 
Cartesian flux tensor, and n is the outward unit normal vector to a cell face. The cell volume 
and cell face area are represented by SI and a, respectively. The quantity h • n is the flux vector 
normal to a cell face. The tensor h can be written in terms of its components in the x, y, and 
z directions as 


where 


h - (E; - Ey) x + (F/ — Fy')y + (G/ - Gy)z 

(2.2) 

q = [p, pu, pv, pw, E] T 

(2.3a) 

E/ — [pu, pu 2 +p, puv, puw, (E+p)u\ T 

(2.3b) 


3 



F/ = [pv, puv , pv 2 +p, pvw, ( E+p)v] T (2.3c) 

G / = [pw, puw, pvw, pw 2 +p, (E+p)w] T (2.3d) 

Et! = [0,T xx ,T x y,T xz ,UT xx -i-VT X y-{-WT xz — q x ]^ (2.3e) 

Ft, = [0 ,T X y,Tyy,Ty Z ,UT X y + VTyy + WTy Z —q y ]T (2.3f) 

Gt> = [0( T XZ , Tyz, T zz , UT XZ + VTy Z -\~WT zz — Qz\^ (2.3g) 


All quantities have been nondimensionalized as shown in the following equations: 
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The value L is a user-defined reference length. The equation set is closed by using the ideal 
gas law. 

Expressing equation (2.1) in iinite volume form for a single six-sided cell in the domain gives 


( n ) t , k ( H »- l p.j.k °i- l /2,J.k H i+ I hj,k 


*+ 


+ H* 


* J- { hM %k a i,j+ l /2,k + H 1 j,k- 1/2 °Uj,k- l /i H i' j,fc+ 1/2 a \ j,k+ (2.4a) 


A shorthand index notation that is used enables equation (2.4a) to be expressed 


as 


r f 'St 

A,) ' • = ( n 


lit T. [(HVV,*-(H Hfl4 ] 

J p=i,j,k 


(2.4b) 


where i5q - (q 7,+1 - q n ) is the change in q per time step, H is the flux vector normal to a cell 
face, St is the time step, and n is the time level. The lower case subscripts in equations (2.4) 
represent cell-centered values, unless offset by one half, then they represent values at the center 
of a cell face. A diagram of the indexing is shown in figure 1. The equations used to compute 
the cell volume, the cell face area, the time step, and the metrics are found in reference 6. The 
asterisk represents the time level and may be at either the n or //. -f 1 level, depending on the 
cell face being evaluated (ref. 6). Since H is comprised of inviscid and viscous contributions, it 
can be expressed as 


H* — Hy(q L , q R ) + Hy(q L ,q /J ) 

The first-order inviscid flux is computed with Roe’s flux-difference splitting to derive its 
upwind nature. The flux difference at a cell face comes from the solution of a Riemann initial 
value problem posed by Roe. The initial data for the problem are found at the two cell 
centers to which the face is common. The flux difference is split into two parts, one associated 
with the positive eigenvalues of the unsteady inviscid Jacobian matrix and the other with the 
negative eigenvalues. The speed and direction of information propagating toward a cell face are 
directly proportional to the absolute magnitude of the eigenvalues and sign of the eigenvalues, 
respectively. This is to say negative eigenvalues send information from the right toward a cell 
face, whereas positive ones send information from the left, and how quickly this information 
travels is related to the eigenvalue’s absolute magnitude. Roe’s method correctly interprets 
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this wave motion relative to a cell face and computes a flux which correctly models the physics 
of the flow. The flux for a supersonic flow with all positive eigenvalues is solely comprised 
of information from the left, and a supersonic flow with all negative eigenvalues is comprised 
solely of information from the right. If the flow has mixed eigenvalues, the flux is comprised of 
information from the right and left. 

In general, an inviscid flux at a cell face can be written as 

H}(q L , q*) = \ [aH}(qL) + bHj( q*)] - ^D*(q L , q*) (2.5) 

where L and R denote the indices of the cell centers to the left and right, respectively, of the 
face being evaluated. The functions, a and b, are flux weighting parameters defined in terms 
of cell volumes as 


_ 2Qg 20 ^ 

a n R + si L nR + nL 

The parameters, a and 6, lessen the effects of grid stretching, and their formulation is empirical. 
The variable D*, shown in equation (2.5), is the previously mentioned flux difference. It can 
be written as 


D*(q L ,qj?) = M |A| M -1 Aq* = |A| Aq* (2.6) 

where A is the flux vector Jacobian. In equation (2.6) the Jacobian is also shown in terms of 
its right (M), and left (M -1 ) eigenvectors, and a diagonal matrix of its eigenvalues (A). The 
matrices A, A, M, and M -1 are given in appendix A. 

With equation (2.6), the first-order inviscid flux at a cell face is 


1 


H/(qi,qft) = » [oH}(qi) +6Hj(qjl) - |A| (qR,q L ) Aq 


(2.7) 


The inviscid flux shown in equation (2.7) can be thought of as being composed of a second-order 
approximation to the flux at a cell face (first two terms) minus a dissipation term (remaining 
term). If this dissipation is not included, the algorithm is equivalent to a centrally differenced 

algorithm. . . 

For second-order accurate solutions with Yee’s STVD approach, D in equation (2.5) is as 

follows: 


D* 


(q L , qi? ) = M 2 |A| 2 \m 2 1 Aq 2 - minmod(M 1 1 Aqi, M 2 1 Aq 2? M 3 1 Aq 3 ) 


(2.8) 


The subscript 2 references the face at which the flux is being computed; 1, the face behind, and 
3 the face ahead. The min mod function compares differences in characteristic variables at 
these locations and chooses the smallest in absolute magnitude if the signs of the values are the 
same or zero if the signs are different. The first term in equation (2.8) is the first-order term, 
the second is a correction which makes D* second-order accurate if the signs of the differences 
in characteristic variables are the same. 

From reference 6, the viscous stress in a direction normal to a cell wall can be expressed in 
terms of its x, y, and z components as 


^nd 


Re 


+ A 

Re 


( U£ ■ V£) + ( Ur) • Vtj) + (u C • VC)] 

U£d + Ur,Vd + (VC • n) + dr, (Vr/ • n) + d C (V( ■ ft) 


(2.9) 
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where d — (x,y, 2 ), and d is the Cartesian velocity component ( d — u , t?, or w) corresponding 
to d . If the thin-layer assumption and Stokes hypothesis are invoked, the thin-layer viscous 
stresses on a cell wall with a unit normal n become 

r nd = Re \ x + ^ x ' (2.10) 

where \ represents the computational coordinates (£, ?/, or () ) . For the thin-layer approximation, 
X is evaluated in the body normal (<) direction only. Including the remaining directions (£ and 
V) is equivalent to neglecting only the cross derivative terms from the stress tensor. 

From reference 6, the heat flux through a cell wall can be expressed as 


/ry 

Pr 


5£(V£ • n) + e r/ (V ?7 • n) + e^(VC • n)j 


(2.11) 


Again, if the thin-layer assumption is made for the heat flux, equation (2.11) becomes 


* = ]5" [ e x(^X • »)] (2-12) 

where the treatment of \ is identical to that of the stress terms. The viscous portion of H is 
shown as follows: 


Qi?) ^0) ^nx > ^ny 5 T n 4 S ~b (2.13) 

If variables with the superscript * contain information referenced to the i, j< k cell center, they 
are linearized by the following equation: 


K* = K n + 





The variable K is a dummy variable which represents the value being linearized. 

Upon substituting equations (2.5) and (2.13) into equation (2.4b) and performing the 
requisite linearizations as described in reference 6, the governing equations can be written 
as 


/ + 


bt 

20 


£ 

P =i,j,k 


A| + 2- 


m v\ (, , 9HJ \ 

* <Tp+1/2+ ( |A|_2 ^r) 


p- 


M 


>,j,k 


= (^j £ ({[aH"(q/J -(- feH;'(q/;) - D n (q t ,q*)] + 2H" (q t ,q w )} p 

\ / i.j.k p= ,j x k 


'k a p~ l /2 


{[ a H"(qi) -)- bH’l (q/j) - D”(q£, q/j)] + 2H^(q/,,q fl )} p+1 ^(r p+1 ^j (2.14) 

where the inviscid and viscous Jacobians are given in appendix A. 


2.1. Boundary and Initial Conditions 

Pseudo (imaginary) cells are used to specify all boundary conditions. The location of these 
cell centers is 1/2 cell beyond the boundary. The volume of each pseudo cell is set equal to the 
volume directly across the boundary being considered. 

At the body pseudo cell centers, Cartesian velocity components are specified with no slip 
conditions, the wall temperature is given, and a normal pressure gradient equal to zero is 
employed. Because the boundary layer is resolved to have a cell Reynolds number (eq. (4.1)) 
of order 1, the distance of 1/2 cell is very small, and as a result, the specification of the wall 
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boundary conditions at these locations appears to have a negligible effect on the ^solution (ref. .6)- 
Values for the free-stream pseudo cells are taken to be the free- stream quanti .e . Reflect ve 
boundary conditions are used in the upper and lower symmetry plane pseudo cells and in the 
axis-singularity pseudo cells. Pseudo cell values at the supersonic outflow boundary are s 
equal to the computed values across the boundary. For this specification of pseudo cell values 
to be accurate, the streamwise flow at the outflow boundary must be supersonic. Although 
flow in the boundary layer is not totally supersonic, the use of a supersonic outflow boundary 
condition does not appear to adversely affect the solution in this regiom Additional information 
on the boundary conditions for the globally relaxed algorithm is in reference 6 
The flow field is initialized with the free-stream values occupying al the cell 
the computation starts, it is as if the body instantaneously materialized in the free-stream 

flow. 


2.2. Space Marching Algorithm 

The original code globally sweeps in planes moving from the body toward the shock and 
vice versa updating the dependent variables at the cell centers of the plane before moving 
to the next plane. This scheme can be thought of as being point Jacobi within the plane and 
Gauss Siedel outside the plane, since values used from within the plane are at the previous dum 
step whereas values used from outside the plane are the most recently computed. To enable 
marching in the streamwise direction, the global sweep direction was changed to thestreamwise 
direction and the capability of performing multiple iterations (local iterations) within a plane 
before preceding to the next was added. In this state, the code can be used as a space marchi g 

code provided 

1. A starting solution is used to begin the computation 

2. Each marching plane is converged by using local iterations before moving to tin next 

3. Only one global pass is used 

Solutions from this initially modified code were obtained for a 15° sphere-cone and can be found 

in section 4. This modification is referred to as “MODI.’ 

From MODI, efforts were made to reduce the computer resource requirements and make 
the code resemble a PNS solver. The number of axial stations in the streamwise direction was 
reduced from the number required to define the entire geometry to five, the number needed for 
space marching Grid generation, which facilitated the space marching solution procedure was 
induded and is addressed in section 3. The dissipation portion of the inviscid flux was i modified 
to enable the computation of second-order accurate fluxes in the streamwise direction. T 
issues as well as boundary and initial conditions and the handling of the streamwise pressure 
gradient, are discussed in the following sections. The code incorporating all modifications has 
been tagged “LAURA-SM”; the additional SM denotes Space Marching. 

£ 2 1 Boundary conditions. The space marching algorithm also employs pseudo cells 
to specify boundary conditions. The methods used to determine pseudo cell values for the 
upper and lower symmetry planes, free stream, and the body are identical to those used in the 
globally relaxed algorithm. The axis singularity boundary condition is replaced with an inflow 
boundary condition. The values at the inflow boundary occupy actual rather than pseud 
cells At the start of the space marching procedure, the inflow boundary is part of the starting 
solution As the solution proceeds down the body, inflow boundary information is saved and 
carried along until no longer required. When updating dependent variables a the P station wrth 
first-order accurate streamwise fluxes, only inflow boundary information at thep- 1 statmn^s 
required but for second-order accurate fluxes, information at both the p - 1 and p 2 statio 
is required The stations are shown in figure 2. The values at the p + 1 station occupy a 
position equivalent to that of the pseudo cells for the outflow boundary in the globally relaxed 
algorffhm However, in the space marching algorithm, actual cells and not pseudo cells are used 
fo? values at the p + 1 station. Because the location p + 1 is treated like a supersonic outflow 
boundary, the values at the p + 1 station are specified with the outflow boundary condition use 
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in the globally relaxed algorithm. As before with the globally relaxed algorithm, this boundary 
condition is valid only for streamwise supersonic flow. In section 2.3.2, it is shown that, when in 
the subsonic portion of the boundary layer, this specification of the outflow boundary condition 
causes the pressure difference across a cell to be a fraction of the total difference and enables 
stable space marching. 

Each data plane is separated by a body normal grid system. Initially, grid coordinates at 
the p - 72, p - 3 / 2 , p - 'fa, and p + V2 stations are provided with the starting solution, and 
the p + 72 face is computed within the space marching code prior to beginning the solution 
procedure. Just as with the information needed for the inflow boundary, the required grid lines 
are saved and carried along as the solution proceeds down the body until they are no longer 
needed. The starting solution is computed by the global relaxation algorithm (LAURA). 

2.2.2. Second-order fluxes. To compute second-order accurate fluxes using Yee’s STVD 
differencing, differences in characteristic variables (M -1 Aq) across a cell face must be defined 
not only at the face at which the flux is desired but also at the face ahead of it and the 
face behind it. All information necessary to compute second-order accurate fluxes in the 
circumferential and normal directions is available. However, in the streamwise direction, a 
problem arises because of the use of the supersonic outflow boundary condition. Referring' to 
figure 2, the use of this boundary condition results in the dependent variables at the p + 1 and 
p locations being identical, and hence the difference in characteristic variables across the p+ >/ 2 
face is zero, since Aq across this face is zero. If the difference in characteristic variables is zero 
at the face where the flux is desired, the face ahead, or the face behind, the min mod function 
returns a value of zero. Because the difference in characteristic variables across the p + i/ 2 face 
is zero, and this value is input to the min mod function when computing a second-order flux 
at both the p + >/ 2 and p - >/ 2 faces, the value returned by the min mod function is zero. This 
results in the streamwise flux being first-order accurate. This is true for the MODI solutions. 

To avoid this condition, differences in the characteristic variables across the p + i/ 2 face 
must not be used in the min mod function. To compute a second-order flux at the p — i/ 2 face, 
the mm mod function compares differences from the p - 3 / 2 and p - i/ 2 faces only. The value 
returned by the min mod function for the flux at the p - >/ 2 face was also assumed to be the 
value returned by the min mod function for the flux at the p + i/ 2 face. Although this method 
ensures that the streamwise flux at a cell face will not automatically be first-order accurate it 
no longer assures that the streamwise flux is TVD. 

The formulation of the equations for the space marching algorithm is identical to the global 
formulation. Dependent variables in streamwise planes are updated by using equation (2.18) 
with the above modifications in the min mod function until the solution for the plane is 

converged. This single formulation for the two algorithms enables easy movement from one 
to the other. 

2*2.3 . Streamwise pressure gradient The handling of departure solutions is similar 
in nature to that employed by Vigneron et al. Vigneron multiplies the streamwise pressure 
gradient by a coefficient lj (ref. 16): 


7 

U ~ 1 + (7 - 1 ) M ? 

The value of u>, which enables stable space marching through the subsonic boundary layer, was 
determined through an eigenvalue analysis of the steady inviscid and viscous flux Jacobian. This 
value limits the streamwise pressure gradient so that the inviscid eigenvalues remain positive 
and the viscous eigenvalues remain real and positive. The physical interpretation of u is that 
it represents the fraction of the pressure gradient associated with upstream information. The 
remaining portion (I — u>) is associated with downstream information. In single sweep space 
marching PNS codes, the (1 -u) portion is usually neglected. Including this portion creates an 
ill- posed problem when space marching. In LAURA-SM, the pressure at a cell face in the flux 
vector is monitored rather than the pressure gradient. The pressure at a cell face can be divided 
into contributions associated with downstream and upstream information. To determine the 
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respective contributions, the flux differences must be known and summed according to the signs 

of the Jacobian eigenvalues. . .. . 

Within the formulation of an upwind-biased flux computed with flux-difference splitting 
is the knowledge that the flux contains portions of upstream and downstream information. 
The value k which represents the contributions to the pressure at a cell face from upstream 
information, and the remaining portion (1 -k), which represents information from downstream, 
occur naturally within the formulation of the upwind-biased momentum flux at a cell lace unlike 
Vigneron’s t j, which must be computed. 

As an illustration, consider the first momentum-flux term in the streamwise flux vector. 
The characteristic variables (M _1 Aq) are 


n = c 2 (p R - p L ) + {p r - pl) 

(2.15a) 

r 2 = Pr 2 Pl 2 ( V R - V l) 

(2.15b) 

r 3 = pfpfiWR - ^l) 

(2.15c) 

T\ = 2 p[ 2 Wr ~ Ul) + (pr _ Pl) 

(2.15d) 

r '5 = ('Pp 2 p'/. 2 L -Ur) + (Pr _ Pl) 

(2.15e) 


where 


U — un x + vriy 4- wn z 
V — ul x + vly 4- wl z 
W = um x 4- vrriy + wm z 

and c is Roe’s averaged speed of sound as defined in appendix A. The variables with subscripts 
L and R represent variables at the cell center to the left and right of the cell face, respectively. 
Variables with the subscript L can be thought of as upstream information, since they are 
upstream of the cell face, whereas those with R can be thought of as downstream information. 
Multiplying through by the absolute value of the eigenvalue matrix and the second row ol t e 
right eigenvector matrix M and isolating the pressure contributions gives the pressure difference 
associated with the second element of the vector D*(q^,qz,) in equation (2.6) as 

A p£ = [M x { 1 - M) + Mn x ](p R - p L ) ( 216 ) 


where 

M = - 


u 

M x = - 
c 

(U and u are computed with Roe’s variables.) The derivation of equation (2.16) is given 
in appendix B. Equation (2.16) was formulated by assuming subsonic flow in the streamwise 
direction and is only valid for Mach numbers from 0 to 1. For Mach numbers outside this range, 
equation (2.16) must be recomputed to account for the different signs of the eigenvalues. 

The pressure at a cell face is 

p(qi,>qb?) = 9 ipr + 0 2Pl ( 2 ' 17 ) 
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where 


ff l = \ K 1 - M)n x + ( M - 1 )M X \ (2.18a) 

°2 = ^ [(M + 1 )n x + (1 - M)M X ] (2.18b) 

Replacing n r with n y or n z will yield similar expressions for the remaining two momentum 
fluxes. 

For one-dimensional flows, equations (2.18) become 

= l -(M - l) 2 = 1 - « 

^2 = 1- ~(M l) 2 = K 

The pressure at a cell face can then be written in terms of k as 

P( c Ili *r) = (! ~ *)Pr + *Pl (2-19) 

where k represents the fraction of the upstream influence and 1 - k represents the fraction of 
downstream influence. 

Examining k as a function of Mach number shows that it is fully consistent with the upwind 
formulation, and as such, the correct physics is reflected in the pressure. For a Mach number 
equal to 0, the pressure at a cell iace is composed of equal fractions of upstream and downstream 
pressure, whereas for a Mach number of 1, the pressure contains no downstream influence. 

The streamwise pressure difference across a cell is 

dp 

— = - [(1 - k) Pi + K Pi -l}i- 1/2 + [(! - «)p <+ 1 + *ft] i+1/2 (2.20a) 

Because of the specification used for the supersonic outflow boundary condition, is equal 
to pi, and enables the gradient to be expressed as 


^7 = K i-l/2(Pi - Pi- 1) (2.20b) 

Although k and Vigneron’s u > serve the same purpose, differences do exist. The value k modifies 
the pressure, whereas u> does the same for the pressure gradient. Because of this difference, 
equation (2.20b) is not interpreted as the allowable portion of the pressure gradient, but as 
an allowable prescription for averaging the pressure at a cell face. Computations with this 
formulation have demonstrated excellent agreement with the global algorithm and experimental 
data as is shown in section 4. 

3. Grid Generation 

To generate the grid at each p + 3/2 station (fig. 2), the axial step size A z must be 
determined. Within the starting solution, an initial axial step size is given. For the sphere- 
cone, subsequent step sizes are computed by multiplying the previous value of Az by a constant 
factor. 

For the HALLS shuttle shown in figure 3, rather than computing a step size based on an 
arbitrary parameter, it was based on flow variables. In one dimension, the Courant number 
can be expressed as 
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cSt 

Cn= a7 


(3.1) 


where c is the speed of sound, St is the time step, and A z is the step size. Solving equation (3.1) 
for A z gives 

Az = — (3-2) 


C N 


The time step formulation used in LAURA and LAURA-SM is given in appendix D of refer- 
ence 6 as 


St 


K 


(3.3) 


where K is the sum of the absolute values of the inverse transit time for a convective wave to 
cross a cell in the f, tj, and < directions plus the inverse transit time for an acoustic wave to 
cross a characteristic length of a cell. Substituting equation (3.3) into equation (3.2) gives 


(3.4) 


All values used to compute Az were taken from the leeward symmetry plane at the p station. 

The circumferential locations of body points depend upon the angle a (fig. 4). The angle 
a has a value of — n/2 on the windward symmetry plane and 7r/2 on the leeward symmetry 
plane. The angle a, 


er = ^ — 7r( 1 — 0) (3-5) 

through the specification of 0 in equation (3.6), 

<; i>(6 ) = Cl 0 5 + C2# 4 + C3# 3 + C4 0 2 + c§9 + / (3-6) 

controls the circumferential spacing and clustering of body points. The function 0 itself specifies 
a percentage of points to be distributed up to the corresponding value of 8 , and 0 (< dcp/d6 ) is 
related to the degree of clustering at that location. Values of 0 that give the desired values of 
the angle a are determined from the constraints used to solve for the constants (ci, C2, C3, C4, C5) 
in equation (3.6). The constraints were determined a priori through a trial-and-error process. 
The specific constraints for the sphere-cone are 

0( 0) = 0.0 <(>(0.5) = 0.5 <X1) = 1-0 

0'( 0) = 1.0 <t>'{ 0.5) = 1.0 = 1.0 

As 6 varies linearly from 0 on the leeward surface to 1 on the windward surface, these conditions 
give uniform spacing around the body. The corresponding constraints used for the shuttle-like 
geometry are 

0(0) = 0.0 4>(8) = 0.5 0(1) = 1.0 

0'(O) = 4.0 0'(0) = 0.4 0'(1) = 1.5 

The value of 8 in these constraints corresponds to the location of the orbiter’s wingtip, which for 
each station must be determined through an iterative process. These constraints cluster points 
near the windward symmetry plane and wingtips. A variation of the true shuttle geometry was 
used in this work. This variation, which models the windward side of the shuttle, is referred to 
as the HALIS shuttle and is based on the QUICK geometry package (ref. 23). 

The distance of the outer grid boundary (A3) above the body surface for the sphere-cone 
was determined by modifying the exponential stretching function found in reference 17. 

A 3 = A noS e(C 2 C^ 1 + 1.0) (3-7) 
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The equations for C\ and ex are 


'(r+ 1 . 0 ) - 

/ r -f 1.0 ) 
V r - 1.0 J 

ex n 


(r + 1.0) 
\r — l.flj 

ex 

1 + 1.0 



ex =l. 0 -^dy 
2 max 

where z max is the axial length of the sphere-cone and ^body * s the axial body coordinate value 
at the p + 3/2 station. The value of C 2 is used as a growth factor and is the ratio of the shock 
standoff distance at the tail of the cone (A^ii) to that at the nose (A nose ). The parameters r 
and exi are used to control how rapidly A3 grows as the algorithm marches downstream and 
have values of 25.0 and 0.3, respectively. 

To determine the outer boundary for the shuttle grid, a subroutine from the global algorithm 
(ALGNSHK) was modified for use in the space marching algorithm. In the global method, 
ALGNSHK adjusts the grid to the captured bow shock and clusters points near the body. In 
the space marching algorithm, ALGNSHK was modified to estimate the shock standoff distance 
at the p— 1/2 and p - h 1/2 stations. A distance for the outer grid boundary at the p-f 3/2 station 
was then extrapolated from this information. 

A body normal grid system is used. A unit vector normal to the surface is computed, 
and each component of the normal is multiplied by a constant. The constant increases as a 
function of normal distance away from the body such that the points at the p-f 3/2 station are 
distributed proportionally to the points at the p -h 1/2 station. 

4. Results and Discussion 

4.1. MODI— Cone 

The effects of the changes made in MODI were examined by computing the flow over a 15° 
sphere-cone at an angle of attack of 0°. The cone had a nose radius equal to 0.009 meter and 
a length of 70 nose radii. Experimentally, this geometry was studied by Cleary (ref. 26). The 
nominal free-stream conditions are given in the following table: 


Mach number 10.6 

Reynolds number per meter 120000 

Pressure, N/m 2 132.2 

Wall temperature, K 294.0 

Temperature, K 473.0 

Velocity, m/s 1461.0 

Density, kg/m 3 0.00972 


The predicted values of surface pressure and heating from the MODI solutions were 
compared with their corresponding values from the globally iterated solution. In addition, the 
MODI surface heating values were compared with the experimental data of Cleary (ref. 26). 
In all cases, the heating was normalized by a theoretical value of the stagnation point heating, 
and the pressure plotted was the nondimensional pressure. Both pressure and heating were 
plotted against z'/R n, where z' was the axial coordinate value in meters and R\ was the nose 
radius, also, in meters. 

Two MODI cases are presented, one encompassing the full length of the cone (70 nose 
radii) and one covering just beyond the overexpansion recompression region (24 nose radii). 
Both these cases use a body normal grid consisting of 21 equally spaced circumferential points, 
65 points between the body and the free stream, and 27 points along the body. The resolution 
of the boundary layer necessary for accurate heating predictions was achieved by clustering 
points near the body to maintain a cell Reynolds number of order 1. The cell Reynolds number 
is defined as 
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(Re) wall 


/ pc\\ 

U/ 


(4.1) 


wall 


where 1 is the cell height. 

The grid generation routine has not been modified, and as a consequence, the grid used tor 
the first MODI computation has spacing in all computational directions equal to that used to 
obtain the global solution. Attempting to lessen the effect of grid-induced differences between 
the global and MODI solutions, the grid over the nose region from the global solution was used 
as the grid for the starting solution. In doing this, the combination of the starting solution gri 
and the MODI grid is equivalent to the grid used to compute the global solution. 

Solutions were converged with a large Courant number (~100) at the start of the iteration 
process to damp the low frequency errors. After a set number of iterations, the Courant number 
was lowered to damp the high frequency errors. Convergence was also enhanced by the use 
of a spatially varying time step tied to a constant Courant number Solutions are considere 
converged when local iteration has reduced the error of the right-hand-side resudual of e 
governing equations to 1.0E-05. The error is measured by using a L 2 error norm. To reach this 
criterion for the MODI solutions, 600 local iterations evenly split between Courant numbers ot 

100.0 and 3.0 were run. . , . ... . . 

The first MODI solution, covering the entire body, compares well with the global solution. 

The streamwise surface pressure is in good agreement even with the large marching steps 
taken. Figure 5 indicates that pressure agreement is within a fraction of a percent excep in 
the recompression area. In this area, MODI underpredicts the minimum pressure by 4 percent 
compared with the global solution. As the flow recompresses, a 13-percent difference in surface 
pressure is noted in figure 5. Once the streamwise pressure gradient becomes small, the two 
methods are in excellent agreement. 

The heating comparison is shown in figure 6. Figure 6 shows good agreement, excep 
in the recompression region where MODI underpredicts the global heating by 10 percent. 
Agreement with the experimental data of Cleary (ref. 26) is slightly better m figure 7 wit, i an 
underprediction of 7 percent in this area. Over the remaining portion of the cone, MODI and 
the globally relaxed solution compare equally well with the data of Cleary. 

The differences in the recompression region may be due to the marching step being too 
large. As a result, the extrapolated value of k at the p + 1/2 station may have been a 
poor estimate of the amount of pressure to be retained. Also, the difference m the order 
of accuracy of the streamwise fluxes between the space marching (first-order) and the global 
(second-order) algorithms may have made a small contribution to the differences seen in the 
recompression region. Decreasing the step size should help lessen this discrepancy by providing 
better resolution in the overexpansion recompression region. This is examined in the second 

MODI case. . , . , , • c „i.;„i 

The second MODI solution was computed by using a grid with the same circumferent 

and normal spacing and one third the streamwise spacing as the first. The streamwise spacing 
of the first MODI grid varied from 0.004 to 0.058 meter. The second MODI grid varied from 
0.001 to 0.019 meter. For both grids, the spacing increased exponentially. The grid dimensions 
and convergence method were unchanged. 

Solutions for surface heating and pressure are much improved over the first MODI predic- 
tions. The difference in the minimum pressure is only 2 percent, and figure 8 also indmates 
that the increased resolution of the overexpansion recompression region has removed the s l 
in pressure experienced in the first solution. Figure 9 reveals improved heating predictions 
when compared with the global values. Unlike the first MODI solution, good agreement is 

obtained in the overexpansion recompression region. 

In comparing the two MODI solutions, the overexpansion recompression region is apparently 
more sensitive to step size than other parts of the flow. Poor agreement in this area is found 
in the first MODI solution, and good agreement is obtained in the second solution with the 
only difference being the step size. Outside of this area both solutions compare equally as well 
independent of step size. 
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4.2. LAURA-SM— Cone 


The flow over the cone described in section 4.1 was recomputed with LAURA-SM to examine 
the effects of the additional modifications. Again, surface pressure and heating are compared. 

Four cases were run. Each had a different factor by which the axial step size grew linearly. 
Factors were varied to see the effect of the step size on stability and solution quality, and 
to determine in terms of computational time whether a large or small marching step is more 
efficient. The initial step size for all runs was 0.0012 meter and was chosen based on the second 
MODI case. The four axial step size stretching factors (FKT) examined were 1.006, 1.10, 1.14, 
and 1.17, corresponding to 233, 45, 35, and 31 solution planes. Factors beyond 1.20 were not 
explored, since this would lead to excessive grid stretching in the streamwise direction. 

The effect of step size on surface pressure is shown in figure 10. In general, excellent 
agreement with the global solution is found, except at the end of the recompression region. The 
expanded scale on the figures indicates that, as the stretching factor is increased, a deviation 
between the global and space marching solutions at the end of the recompression becomes 
prominent. The larger the stretching factor, the longer the delay in recognizing the end of the 
recompression. A maximum difference of 2 percent is associated with the largest factor (FKT = 
1.17) and a difference of less than a percent with the smallest factor (FKT = 1.006). In looking 
at the location of the minimum surface pressure, decreasing the stretching factor produced a 
lower minimum pressure and moved its location downstream. Comparing the solution computed 
with the largest factor with that of the smallest shows that the minimum pressure is lower by 
1/2 percent and has moved downstream approximately 1 nose radius. 

The effect of step size on surface heating and comparisons with the global solution are 
shown in figure 11. Overall the space marching solutions compare well with the data. The 
predictions for surface heating from the space marching code more closely match the predictions 
from the global relaxation algorithm than the experimental data of Cleary (ref. 26). As the 
stretching factor is decreased, differences between the global and space marching heating values 
are resolved to less than 8 percent. The space marching solution also shows a small inflection in 
the heating just downstream of the minimum in pressure. The inflection becomes less noticeable 
as the step size is reduced. As did the location of the minimum in pressure, the inflection in 
heating also moves downstream as the stretching factor is decreased. Both the changes in 
location and magnitude are small. In comparing results from the largest and smallest factors, 
a 1-percent difference in the minimum value exists, and the location moves downstream 1 nose 
radius. 

The deviation noted between the space marching solutions appears to be step size related. 
Increased resolution of the flow in the streamwise direction resolves these discrepancies but 
adds to the computational time to obtain a solution. Because the global and space marching 
solutions are in good agreement outside the recompression region, regardless of the step size, 
including the overexpansion-recompression region in the starting solution for flows where it 
occurs early and is well defined may enhance computational efficiency, since larger step sizes 
can be used outside of this region. While a larger step requires more iterations to converge 
than a smaller one, the pass over the vehicle is performed in less computer time with the larger 
step, as is shown in the next paragraph. 

Two advantages of using a space marching algorithm are that less computer time is needed 
to obtain a solution when compared with a globally iterated algorithm and the space marching 
algorithm can be used on a wider range of computers. The following table shows the computer 
time required to compute the global and space marching solutions: 


Code 

Time, sec 

Factor 

Steps 

LAURA-Global 

12 800 

1.00 

39 

LAURA-SM (FKT = 1.006) 

7429 

1.72 

233 

LAURA-SM (FKT = 1.100) 

7880 

1.62 

45 

LAURA-SM (FKT = 1.140) 

7256 

1.76 

35 

LAURA-SM (FKT = 1.170) 

6949 

1.84 

31 
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The time shown is the CPU time required for the VPS-32 at the Langley Research Center The 
first coTumn fists the code used, the second the CPU time, the third fists the ratio of the global 
to space marching CPU time, and the fourth shows the number of streamwise steps taken. The 
CPU Times for tta space marching solutions include the 3102 -c required .u » obtarn the start ng 
solution The starting solution accounts for approximately 45 percent of the CPU tune in 
space marching cases. This is due to the subsonic region at the nose. It is possible to red 
the' computational time of the space marching and global algorithms approximately 40 percent 
by freezing the left-hand side of the governing equations, and updating it periodically (ref. 6). 

4.3. LAURA-SM — Shuttle Orbiter 

To evaluate how the space marching algorithm performs when exposed to thvee- 
dimensional effects the flow over the HALIS shuttle geometry was computed. This geometry 
models the windward surface of the Space Shuttle orbiter. The HALIS model fills in the leewar 
surface with elliptic cross sections. Figure 12 shows both geometries. 

Predicted surface heating values from the space marching algorithm arc ' 
unpublished experimental heating data of Miller. At higher angles of attack. Miller s data have 
been published (ref. 27), but are not applicable for use in this instance since the wmdw 
surfece^s subsonic . The flow was computed to correspond to the wind tunnel conditions used 
by Miller. The nominal flow conditions are as follows: 


. 10.0 
690000 
110.7 

300.0 

492.0 
1410.0 

0.00784 
. . 10 


Mach number 

Reynolds number per meter 

Pressure, N/m 2 

Wall temperature, K 

Temperature, K 

Velocity, m/s 

Density, kg/m 3 

Angle of attack, deg 

It should be noted that in Miller’s experiment, the actual shuttle geometry was used and 
not the HALIS shuttle geometry shown in figure 12(b). Even though the leeward surface 
used to obtain the predicted results is different for that used in the experiment, comparisons of 
windward surface values can be made. The circumferential flow around the dune. land r«™i£P» 
is supersonic, and minimal influence is transmitted circumferentially through the boundar> 
layer therefore, the leeward surface shape has negligible effect on windward values. 

The space marching solution begins at z'/z ie{ = 0.062 and continues toward a value 
of 0.942. The axial distance of the experimental model from the nose to the base z ref is 
0 196 meter. The space marching algorithm computed values at 97 streamwise s a lons _ 
Courant number of 50.0 was used to converge each station. A station was considered converged 
when locaHteration had reduced the L 2 norm below a value of 0.2E-06 A change in the grid 
system was also made at z'/z tei = 0.350. At this location, the grid system was changed from 
body normal to axis normal. The change was made to circumvent the occurrence of normal 

erid lines crossing in the area where the wings begin to flare. 

The heating predictions agree within 5 percent when compared with the experimental data 
of Miller (ref. 27). Figure 13 shows the predicted values following experiment, until predicted 
values begin to diverge at a location z/z ref of 0.55. At this location, a cross-flow Ao^ontibe 
leeward surface intensifies. Beyond this location, the heating agreement worsens. The space 
marching algorithm fails at 0.70, where the algorithm encounters reversed streamwise flow o 

the leeward surface near the symmetry plane. , . . 

Two potential causes of this divergence have been identified^ First the developing cross 
flow shock is oblique to the grid. Such orientations cause enhanced dissipation across the 
captured shock and may induce physically incorrect pressure fields in the ‘^mediate « 
that can induce separation. Limitations in the present grid generation capabilities P^c d 
resolution of this problem. Second, the implementation of the second-order corrections in Yee 
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STVD scheme must be modified in the streamwise marching direction. The min mod function 
used in the STVD approach selects the gradient of characteristic variables which is guaranteed 
not to introduce any new maximums or minimums into the problem or revert to first-order 
accuracy^ By doing this, shocks and other discontinuities are captured without oscillations. 
I he modification made in the streamwise direction, which excludes extrapolated gradients 
at downstream cells from consideration in the min mod function, may degrade the solution. 
This change may neutralize the TVD characteristics of the original scheme, particularly in the 
subsonic portion of the boundary layer. However, this restriction cannot be removed without 
sacrificing second-order accuracy in the marching direction. 

Pressure contours for the global and space marching algorithms at three cross-sectional 
locations are shown in figure 14. The presence of the cross-flow shock is confirmed in these 
contours. The patterns of the contours are similar for the global and space marching solutions. 
Both indicate a high gradient region near the wingtip and the formation of a cross-flow shock 
near the leeward symmetry plane. A thickening of the shock is also noted in the space marching 
contours as the solution proceeds downstream. 

5. Concluding Remarks 

The Langley Aerothermodynamic Upwind Relaxation Algorithm (LAURA) code has been 
converted to perform space marching of the thin-layer Navier-Stokes equations in the spirit of a 
parabolized Navier-Stokes (PNS) solver. This study has demonstrated that the intrinsic upwind 
treatment of the streamwise pressure combined with an appropriate outflow boundary condition 
avoids departure solutions and stable marching is possible without the explicit inclusion of a 
limiting variable within the governing equations. Excellent agreement with the globally iterated 
algorithm for surface pressure was obtained with the new code on a sphere-cone geometry, and 
computed heat-transfer distributions on the sphere-cone and the HALIS shuttle orbiter were 
m good agreement when compared with experimental values. 

Care must be exercised when computing through an overexpansion-recompression region. 
Taking too large a step in this area may not cause the solution to diverge but to yield an 
inaccurate solution. Unfortunately, a prescription for the best combination of initial step size 
and stretching factor which would give the best solution and minimize computer time is not 
known. The best combination must be determined through experience. 

The LAURA-SM algorithm was faster than the globally relaxed algorithm, but only by a 
factor of approximately 1.6. Room for improvement may lie in including regions which require a 
small step size, such as the overexpansion- recompression region associated with the cone, in the 
starting solution and initiating space marching beyond the overexpansion-recompression region. 
Another approach may be to freeze the inverted implicit part (left-hand side) of the governing 
equations and update it periodically. Doing this can reduce computer time by approximately 
40 percent and would have no effect on the converged solution. This same freezing procedure 
can be used with the global algorithm, saving a similar percentage of computer resources 

The space marching HALIS shuttle orbiter solution fails at the cross-flow shock on the 
leeward surface. Since the second-order streamwise flux is no longer total variation diminishing, 
it may introduce spurious oscillations in the solution which cause divergence. If the second- 
order streamwise flux is at fault and the space marching formulation of the governing equations 
is to remain unchanged, a way of specifying the gradients of characteristic variables consistent 
with the global approach must be developed if the min mod is to function as intended. 
However, if second-order accurate streamwise fluxes are to be computed, this is not possible, 
since downstream gradients of characteristic variables are required in subsonic regions and are 
unavailable. The only other way around this condition is to use first-order streamwise fluxes 
and sacrifice overall second-order accuracy. Another factor which may have caused the space 
marching algorithm to fail is grid alignment. Not aligning the grid with a shock can result in 
the pressure field being different from when the grid is aligned. The resulting pressure field 
from the space marching solution could have caused the reversed flow found on the leeward 
surface. This question was not resolved because the required adaptive-grid capabilities are not 
incorporated in this code. 
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Appendix A 
Jacobian Matrices 

The unsteady inviscid flux Jacobian matrix A is 

0 n x n y n z 0 

-uU + /3an x U + un x - un x 0 un y - vn x 0 un x - wn x 0 0n x 

-vU + 0an y vn x - un y 0 U + vn y - vn y 0 vn z - wn y 0 0n y 

-wU + 0an z wn x - un z 0 wn y - vn z 0 U + wn z - wn z 0 0n z 

U{a0-H) Hn x — uU 0 Hn y -vU0 Hn z -wU0 7 U 

The matrix A can be decomposed into it right (M) and left (M" 1 ) eigenvectors and a 
diagonal matrix of its eigenvalues (A). These matrices are 





• a0 - c 2 —0u -0v —0w 0 

-V l x ly l * 0 

M' 1 = -W m x rn y rn z 0 

a0 -Uc cn x — 0u cn y - 0v cn z - 0w 0 

. a0 + Uc —cn x — 0u —cn y — 0v —cn z — 0w 0 


■u 0 0 o o 

QUO 0 0 

A= 0 0 U 0 0 

0 0 0 U + c 0 

. 0 0 0 0 U - c 


(A2) 


(A3) 


(A4) 
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where 


U = un x + vn y + wn z 
V = ul x + vly + wl z 
W — um x + vTriy + wrn z 

a = ~(u 2 + v 2 + w 2 ) 

c 2 = {H- a)0 

In order for equation (2.6) to be an exact equality, specially averaged variables known as 
Roe’s averaged variables are used in the computation of the matrices in equations (Al) through 
(A4). These variables are defined as 


u = ai(a 2 u R + u L ) ' 
v = ai(a 2 v R + v L ) 
w — a\{a 2 w R + wi) 
H = a\(a 2 H R + H L ) , 


where 


a i 


pl ) 


1/2 


a 2 = 


1 4- a\ 


The varriable (3 is defined in terms of the ratio of specific heat as 


(A5) 


0 = 7 - 1 


The viscous Jacobian, shown in equation (A6), was used in equation (2.14) and was taken 
from reference 6: 

dHv 

dq 

where 

*1,1— S = 0 

* 1 — 4,5 = 0 

m = -(w+^) 

*5,1 = «&2,i + vfy'i + wb^i - p- [£■ - (u 2 + v 2 + w 2 )] 


^(Vx • n) 

Re p 


[**»>] 


(A6) 


b 2,2 


n x n x 

3 


b 3,2 = & 2,3 = 


3 
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n x u z 

H,2 - b 2A = ~ 


^ 4,3 = ^ 3,4 = 

&3,3 = 1 + 
^ 4,4 = 1 + 


TlyTlz 

1 “T” 

n v n v 

3 

n z n^ 

3 


7 U 


f>5,2 = ^62,2 + U&3,2 + w ^ 4,2 ~ 


7 v 

65. 3 = ub 2y 3 + u6 3)3 + ^ 4,3 - p^ 

7W 

65.4 = ub 2,4 + ^63,4 + ^64,4 - p^; 


All variables are values computed at cell faces with the exception of u, v, and w, which are 
cell-centered values. 
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Appendix B 

Streamwise Pressure Gradient 

Equation (2.16) is computed with equation (2.6), which is 


D(q x, q/{) = (MA+M- 1 - MA~M _1 )(q/j - q L ) 
where A + and A - are given by 

^ + |A| _ A - |A| 

2 2 


(Bl) 


(B2) 


The first term in equation (Bl) represents information associated with the positive eigenvalues 
and the second term is associated with the negative eigenvalues. For stable space marching, 
reverse flow is not allowed; therefore, the streamwise Mach number must be greater than or 
equal to 0 throughout the flow field. In the boundary layer, the streamwise Mach number is 
further restricted to be less than 1, since the flow is subsonic. Accordingly in the subsonic 
boundary layer, four eigenvalues will be positive ( U,U,U,U + c) and one will be negative 

(D — c). Knowing the eigenvalue signs, equation (Bl) can be expanded accordingly as shown 
below. 


-zv 


(U+\U\)(c 2 Ap + Ap) 

D 2 

= MxI 

(i U + \U\)p l l 2 pf AV 

D 3 

(u+mp^pf aw 

Da 

-D b . 


[(17 + c) 4- (|£7 + c|)j {cp X [ 2 p£ AU + Ap) 



0 


0 

0 

0 

0 

[(U - c) - (| U- c|)] {^cpfpf A U + A p) 


(B3) 


Multiplying the vectors on the right-hand side of equation (B3) by the second row of M gives 

D 2 = ( Uc 2 A p + U A p) + [Upfpf Av) l z + (up^p 1 ^ 2 AW^j m x 

+ [(U + c) (cp l l 2 pf A U + Ap)] (^|^) - [(U - c) (cp X [ 2 pf AU + Ap)] (^gr) (B4) 

Isolating the pressure terms in equation (B4) gives the same pressure contribution to the 
momentum flux associated with the flux difference. The result is 

a a ( —Uu Un x u\ 

Ap { =Ap(- ?r + - r + -j (B5) 


Equation (B5) can be simplified to 


(B6) 
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= Pr [M x ( 1 - M) + Mn x \ - p L [M x { 1 - M) + Mn x ] 



where 


U u 

M — — M x — ~ 

c c 

Variables with a tilde (~) over them are Roe's averaged variables and should be computed 
as shown in appendix A (eqs. (A5)). The delta (A) indicates a difference across a cell face 
computed with cell-centered values, for example 

A p = Pr- PL 
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Figure 4. HALIS orbiter cross section 
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Figure 5. Global versus MODI surface pressures for z'/Rn - 70, p'^ - 0.00972 kg/m , 1461 m/s, and 

Rn = 0.89 cm. 
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Figure 7. Experimental versus MODI surface heating for z'/Rjv = 70, p'^ - 0.00972 kg/m 3 , S'^ - 1461 m/s, 
Rn = 0.84 cm, and qvpf = 0.2156 MW/m 2 . 



Figure 8. Global versus MODI surface pressure for z'/R N = 30, p'^ = 0.00972 kg/m 3 , S' = 1461 m/s, and 
Rm — 0.89 cm. 
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Figure 9. Global versus MODI surface heating for z'/Rn = 30, p ^ — 0.00972 kg/m 3 , S'^ — 1461 m/s, and 
q rP f = 0.2156 MW/m 2 . 
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Figure 13. Shuttle in windward symmetry plane surface heating for p ^ = 0.00782 kg/m 3 , 5^ - 1410 m/s, 
z re { = 32.66 m, a = 10°, and q Te{ = 687 kW/m 2 . 
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